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We develop a statistical model for a confined chain molecule based on a monomer grand 
canonical ensemble. The molecule is subject to an external chemical potential, a backbone 
interaction, and an attractive interaction between all monomers. Using a Gaussian variable 
formalism and a mean field approximation, we analytically derive a minimum principle from 
which we can obtain relevant physical quantities, such as the monomer density, and we 
explore the limit in which the chain is subject to a tight confinement. Through a numerical 
implementation of the minimization process we show how we can obtain density profiles in 
three dimensions for arbitraty potentials, and we test the limits of validity of the theory. 

PACS numbers: 05.20.Gg, 82.35. Lr, 87.15.A- 

I. INTRODUCTION 

The classical theory of fluids in thermal equilibrium is a highly developed discipline. Along with 
specific physically motivated approximations have come tools of more general validity and utility. 
Studies of polymers have of course been extensive, and useful relationships between polymer melts 
and fluids established [1], as well as between single polymers and fluids. This paper represents an 
initial investigation, aimed at understanding which concepts borrowed from bulk fluid studies [2, 3] 
still remain relevant in a polymer setting. Our focus will be on a single polymer chain confined 
by external forces, and to minimize the needed information input, on single homopolymers. We 
will also attend in the main to idealized models in which the polymer is simply a chain of unit 
monomers of a few degrees of freedom, but will indicate how this restriction can be rewardingly 
removed, on the way to a realistic polymer representation. 

We will aim at both analytic simplicity and reasonable suitability for the ultimately necessary 
computational procedures. For the former, we will work in a "monomer grand-ensemble" [4, 5] 
in which the number of monomers per polymer is distributed, but will show that this need not 
be a drawback. For the latter, we will favor minimum principles to be able to better control 
computations. This does restrict the category of systems to be studied e.g. to purely attractive 



pair interactions, and will be generalized in later work, soon to be reported, as such restrictions 
are removed. 



II. THE REFERENCE SYSTEM 
A. Notation 

Let us be a bit more explicit. We have in mind an ordered chain of N equivalent monomers, 
the j th being specified by its degrees of freedom rj. The order is maintained by a symmetric 
next neighbor interaction potential of Boltzmann factor w(ri,ri + ±), depending of course on the 
inverse temperature [3. Any two monomers can also interact via an interaction Boltzmann factor 
e(ri,rj) = exp(—/3(f)(ri,rj)) and, crucially for the applications we have in mind, the polymer is 
constrained by an external potential u(rj). 

It is now convenient to imagine that the homopolymer in question is both in thermal equilibrium 
and in number equilibrium, i.e. that it is the result of monomer addition and absorption from a bath 
of non-interacting monomers. The reaction equilibrium is analogous to that of the grand canonical 
ensemble for a fluid, but with a significant difference. Suppose that is the monomer canonical 
partition function in its center of mass coordinate system, Q$ of the iV-monomer polymer (N > 1, 
to recognize an object as a polymer). Then if the full system contains N' monomers, N of which 
are bound together, in a volume V, the full system partition function will be: 



Qn '~L w=w Qn - 



(i) 



Since the monomers in the absence of a polymer would have a partition function: 



N'\ 

that attributed to the polymer will take on the form: 



- (p) [Co] = .„ lim_ , n ,„ (m w Q& (3) 



N'\ 

N'-*oo,N"/V&xed (VQ^) 1 

= E CoQ%\ where Co = (N'/V)/Q^l (4) 

N>1 

The obvious analogy with a fluid grand partition function (with the weight 1/N\ excised) will be 
very useful indeed. 

The computation of H^(Co), which generates all thermodynamics and expectations in a thermal 
ensemble, is of course too general to be explicitely solvable, except in very special circumstances. 



Let us therefore start with the evaluation of E^ p \(q) for what may be regarded as the backbone of 
the polymer, that in which the arbitrary pair mutual interaction 0(rj,rj) is set equal to zero. In 
this case, we have at once: 

CoQn = J an)w(r 1 ,r 2 )ar 2 )...C(r N . 1 )w(r N ^,r N )arN)dr N (5) 
where CM = Coe _/3u(r) = e^ r) 

or regarding w(r, r') as the kernel of an integral operator w, ((r)S(r — r') as that of the diagonal 
operator £, and the symbol |1) denoting the vector whose components are 1, 

sW[c] = (i|(r 1 -^r 1 ii), (6) 

subject of course to convergence of the series (4). The corresponding "grand potential" is: 

^)[ C ] = -Iln(l|(r 1 - U ;)- 1 |l) (7) 

and as an immediate consequence the monomer density is given by 

n(r) = -/3C(r)^fiW[C] (8) 
= (1KC 1 - w)- l \r)(r\(C l ~ wY^CH^/m- 1 ~ w)- l \l) (9) 

B. Number Distribution 

The absence of the statistical weight \/N\ suggests an unusual distribution of monomer num- 
ber. To verify this, we first observe from (9) that if N denotes the monomer number in a given 
configuration, than 

N=(N) = Jn(r)dr=(l\(C 1 -w)- 1 C 1 (C 1 -w)- 1 \l)/Z^ (10 ) 

= (C 1/2 I(/ - C 1/2 K 1/2 )^IC 1/2 )/(C 1/2 l(/ - e'M 1 ' 2 )- 1 ^ 12 ) (ii) 
= 1/(1 -A ) + O (1/(1 -Ai)) (12) 

where Ao = Xmax{C l ^ 2 wC}^ 2 ) and Ai is the next largest eigenvalue; Ai < Ao in a confined system. 
If N is large, then 

s W(c) = (c 1/2 |(/-C 1/ V 1/2 )- 1 !C 1/2 ) (13) 

is dominated by the "resonance" at Ao, so that if 

(C 1/2 < 1/2 )V% = AoV% (14) 



with normalized ip, then 

H(p) K]~r^<V%!C 1/2 |i) 2 (15) 

There are now two consequences. On the one hand, we have 

(e ie ") = e *<bS/SCo S (p)[£]/ s (p)[£] ( 16 ) 

= Z^[(e ie ]/Z^[(} (17) 

^ (1 _ Ao )/(l - X e ie ) (18) 

so that the distribution of N is given by 

/(AT) = coef e me (l - A )/(l - \ e i9 ) = A^l - A ) (19) 

the very broad geometric distribution. It would appear that the iV-ensemble must give a very poor 
representation of a given N. But on the other hand, we have from (5) 

Qn = Co-^KC^V^T- 1 !^ 2 ) (20) 

= Co-^-^AolC 1 / 2 ) 2 , (21) 

so that if -F/v is the canonical Helmholtz free energy, 

PF N -paP> = - In Q N + In ~^[C] (22) 
7VlnCo-(A r -l)lnA -ln(l-A ). (23) 

It follows that for any (TV-independent) parameter variation 

0(6F N - 6nW) = (^-L- - <5A (24) 

which vanishes if Co in (14) is chosen so that 

N = 1/(1 - A ) (25) 

We conclude that for this choice of Ao, expectations at fixed TV" and in the monomer number 
ensemble are in fact identical to leading order. 

C. Minimum Principle 

With the confidence that the large spread of monomer number does not detract from the 
usefulness of fi^), we proceed. The original form (6) is simplest to use for the construction of an 



associated minimum principle. It depends upon the fact if K is positive semi-definite, and a and ip 
(not to be confused with ip\, which we will not need again until later) are arbitrary, then according 
to the Schwartz inequality, we have 

(a|lf- 1 |a)(Y>|if|Y') = (K-V 2 a |jf-V2 a } (^1/2^1/2^ (26) 
> (K- 1/2 a\K 1/2 i;) 2 = (a|^) 2 , (27) 

so that (a|if _1 |a) > (°#) 2 / < V'l^lV')- Then indeed 

(aU) 2 , , i, , MW) , s 

Max ^ / 7TF /\ = ag ° at^= , | , a, 28 

and in the context of (6), we conclude that 

= Min^ [ln^lC 1 - w\il>) — 21n<l|V>] • (29) 

III. THE MEAN FIELD STRATEGY 

Our task is now to take into account the mutual interaction </>i(rj,rj) which while typically 
fairly short range in space, can be very long range along the polymer chain. We will confine our 
attention to purely attractive (negative definite) potentials, and to avoid confusion will set 

<h.(ri,Tj) = -4>{ri,rj) (30) 

where (f> is positive definite (as a continuous matrix). The assumption of pure attraction distorts 
local properties of the polymer, but permits large scale effects such as globularity to proceed 
unhindered. 

At first, we need not restrict our attention to polymer chains. Quite generally, if Wq N \t n , fi) is 
the suitably weighted Boltzmann factor for the iV-unit configuration in the absence of (f>, but with 
local chemical potential /x(r), the partition function for the interacting system can be written as 

3^,0] = E / ••• / wl ) N \r N ^)e^' t>ir ^ ) dr N (31) 
N J J 

= E /••• / W^ N \r N ^)e-^n{r i ) e ^ i ^{r i ,r j ) dr N ^ 
N J 

= E / • • • / W<> N \r N ,fi)e-% EjvMjffHrMryMr'Wrdr'^N ^ 
N J 

where J2' omits the i = j contribution, 4>£> is the diagonal part of the matrix 0(rj,rj), and 

n( r ) = S(r - n) (34) 



is the "microscopic" density. The device of Kac, Siegert, Hubbard and Stratonovich [6-9] is to 
represent the Gaussian in (31) (in obvious notation) as a functional Laplace transform 



= J e-^i^e-^'^Dv I J ' e-^'^Dv (36) 

Since J2n I ■ • • I Wq N \v n , ^)dr N = S [/x], eq.(31) can thereby be rewritten as 

E\fi, 0] = J Ho[m - - u]e-§ "•*" lw £) V / J e-% v -+- lv Dv (37) 

with a possible interpretation that the interaction —(f) has been replaced by an ensemble average 
over a fluctuating external field v(r), serving as a sort of graviton shutting back and forth between 
units. 

The kernel of (37) is a Boltzmann factor in field space, and so we may define a field average as 
«GM» = J G[v]E \ M -^ D -v]e-^-1'- lv Dv/E\ fi ,<l,] (38) 

A very suggestive consequence of this notation follows from the observation that for the density 

n(r), 

n(r) = ^^f M (39) 
= - / J^Tf^ ~ ~ v^^Dv /£ [/x, 0] J e-fr+- lv Dv, (40) 

or integrating by parts in v-space (assuming the absence of boundary terms) 

n(r) = - J ^"^(rjHolAi - \<t>D - vje'^'^Dv /Efr, cf>} J e'^'^Dv (41) 

= -r\Kvir)% (42) 

identifying —<pn(r) as the "mean field" {(v(r))). It is the approximate computation of this mean 
field that we must attend to, directly or indirectly. 

If the field Boltzmann factor of (38) is sharply peaked about a function v(r), than of course 
(42) becomes simply 

n (r) = _0-^( r ) ; (43) 
and the field v must now satisfy (Q is again the grand potential — ^ lnS) 

^ (n [» - ^ D - v]j + \v ■ r'v = 0, (44) 
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or simply 



n o(r\n ~ ~ v) = -<p 1 v(r), (45) 



the unsurprising result (comparing with (43)) that n(r) is equal to the "bare" density no in the 
presence of v, to within an additional shift of \4>b- But in addition, we now have in leading order 
approximation the very explicit 

fii (//,(/>) = ^o[« - ^<Pd ~ v] + ■ <\T x v (46) 

which, since (44) represents the minimization of the kernel of (39), takes the variational form 

fii(M> 4>) = Min„O [^ - - v] + -v ■ <?!>~V (47) 

Furthermore, consistency is established by noting that by virtue of (47), fii of (46) implies 

n(r) = = -rMr), (48) 

reproducing (43). 

To apply (47) to the single homopolymer under discussion, we need only insert (29), serving as 
fio> into (47). We then have 



[//,<£] = Min„ ; y, 



'-v ■ <fr x v — 21n<l|V> + ln^le-^'-^ - u;| 



(49) 



where /u'(r) = u(r) — Eq. (49) can be simplified by finding and eliminating the field tp, 

and replacing v by — <fin. To do so we observe that 



6Sh 

Sv 

Integrating over the implicit r, with N = J n(r)dr, then 

^,| e -/%'-tt)| 

{i)\e-w- v ) -w\ipy 



= ^1 = <jr\ + e -^'- v H 2 /^\e-^'-^ - u#) (50) 
ov 



N = — mr> ( 51 ) 



reducing (50) to 

n = Ne-M-^il) 2 /{ip\e-^'-^ (52) 
Since (52) is homogeneous of degree in ip, we are free to adopt the normalization 

(i/;\e-W-»)ty) = N (53) 

converting (52) to 

n = e-M-v)^ (54) 
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(leading back to (53)), and consequently replace (49) by the simple 



(3Qi\p,(f>] = Min r , 



• cf>n - 21n(l|n 1 / 2 el^' + ^) + ln^ 1 / 2 ! (/ - e ^' + ^we^' + ^) 

(r>r >} 



appearing as a /i-dependent density functional, reminiscent of "statistical models" [10-12], of poly- 
mers. Eq.(55) is our main result, valid in the mean field level. The role of the interaction —(f) is 
transparent: the external potential is augmented by the mean field :/x — ► // + <M, and the explicit 
energetic component — \n ■ 4>n subtracted out. Eq.(55) is of course a density functional representa- 
tion with the density profile dependence on \x determined by dropping the implicit (f> and making 
H explicit, by 

<5fii[//,n] /5n(r) = 0. (56) 

It is an extension of the more common form in which one is given a free energy F[n] in terms of 
which fi = 5F(n)/5n, from which it follows that 

Sn(fi,n) /Sn = 0, (57) 

where n] = F[n] — // • n. A corresponding format is obtainable here as well, but at the cost of 
increased complexity. 



IV. NUMERICAL IMPLEMENTATION 



To exemplify the usefulness of the mean field minimization strategy, we have developed a nu- 
merical implementation of the minimization principle (55). 

From preliminary simulation studies we had noticed how in the absence of a repulsive core 
interaction between monomers we could not converge to a stationary density. Instead the density, 
and therefore the total particle number, either kept on growing or shrank to zero, depending on 
the value of the external chemical potential. This numerical instability occures as a result of the 
closeness to the resonance. In these simulations we noticed, on the other hand, that despite the 
fact that the density itself kept changing, the global shape of density profile reached a steady state. 
We therefore decided to study the normalized density p(r) = n(r)/N. Once we had formulated the 
problem consistently in terms of p all the numerical instabilites desappeared and we were able to 
determine the corresponding stady state densities and total particle numbers. 
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In terms of p the minimization functional (55) becomes 



ill [n, <f>] = Min p 



t^p ■ cPp - 21n(l|p 1 / 2 ef V + ln^l (/ - e^'+^^f^'+W)) | p i/2) 

(58) 

Once we have determined p(r), the expression to compute the particle number and the true density 
are 

n{r) = (p(rmi-tl(]p(r>)V 2 ) m 
N = {pir) 1 ' 2 ^ - QwtW) 1 ' 2 )- 1 (60) 

The minimization is done using a Metropolis Monte Carlo algorithm at fixed temperature, 
starting from a suitable ansatz for the functional form of p{r). We can freely impose the form of 
the next neighbor potential, and therefore w(r,r'), the external confining potential u(r), and the 
long range interaction cj)(r,r'). The main numerical difficulty arises from the need of evaluating a 
functional involving 6-dimensional integrals nesting 3-dimensional integrals, such as 

J y^ef^W^OKO^ (6 i) 

It turns out it is prohibitive to try to evaluate these terms directly by the use of a simple grid. 
Already when using 10 discretization points simulations are too slow, while the precision is very 
poor. To achieve higher accuracy in the integration, we use a Legendre-Gauss quadrature method 
with either 6 or 8 points in each linear direction. 

As a test system we have chosen a square well external potential with p = —po inside; a 
harmonic n.n. potential, leading to w(r, r') = exp(— Ch(\r — r'| — pi) 2 ); a long-range step potential 
<p = Ci r @(x — P2), where po, Ch, C\ r are energy constants, and p\ and P2 are characteristic lengths. 
In the following we take the system size to be L = 1, Ch = 5, and p\ = 0.2. As ansatz for the 
probability density we take the composition of three Gaussian in the three spatial directions, for a 
total of 6 parameters: B x ,B y ,B z ,Xo,Yo and Zq. 

P (X, y,z)=A e -B*{x-X f e -B v (y-Y f & -B z {z-Z Q f (62) 

At each iteration of the algorithm the normalization of n is verified so that J p(x, y, z)dx dy dz = 1, 
and the value of A is changed accordingly. 

We first look at the reference system, i.e. zero long-range interaction, figure 1. As expected 
the probability density converges to a single Gaussian centered in the middle of the box, with 
Xq ~ Yq ~ Zq ~ L/2 and B x ~ B y ~ B z . When the attractive potential <f) is turned on the 
Gaussian becomes more peaked (figure 2). ). 



FIG. 1: Left: Interaction between neighboring particles along the chain: harmonic potential only. Right: 
Slice of the probability density profile taken at z = L/2. 




Next, we introduced an asymmetry in the external potential, like a gravitational potential along 
z: li = liq — gz. As expected the minimization of the functional converges to a probability density 
off-centered with respect to z and centered in the middle of the box with respect to x and y (figure 
3). 

Even if the theoretical derivation of eq.(55) is valid only for strictly positive definite <ft, in the 
algorithm we can invert the sign of the long range interaction and make it repulsive, as well as 
attractive. We tried to go beyond the known validity of our approximation, taking a repulsive 
long-range interaction. Our method still holds, provided we consider repulsive interactions which 
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FIG. 3: Probability density profile at z — L/2 in the presence of an external graviational field (g = 1, 
Ci = 2) 
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FIG. 4: Left: Interaction between neighboring particles along the chain: harmonic potential and repulsive 
long-range interaction (C/ = —2). Right: probability profiles at y — z — L/2 with (red) and without (blue) 
the repulsive core. 



are not too strong. As expected a repulsive <f> gives a less peaked density compared to the reference 
system (figure 4). 



V. LIMIT OF TIGHT CONFINEMENT 



We have seen, in Sec. II B, that in the absence of non-neighbor interactions, the large N 
monomer distribution is determined by the "resonance state" ip\ satisfying 

(C V V 1/2 )V% =A<W%, (63) 
11 



where Xmax = Ao ~ 1. Explicitely, since 

-i . . ■,.'„/ : , ,,a -J 



(rliC 1 -™) 11) = (r|C 1/2 (/-C 1/2 < 1/2 ) C 1/2 |l) 

' ' lC 1/2 ^Ao)(V'AolC 1/2 |l) (64) 



1- Ao 

we have under these circumstances, from (9) 



n ( r ) ~ ; 1 7 i>\ ( r f, (65) 

1 — An 



or since ip\ (r) ~ ipi(r), 



n (r) „ pr¥_ = N^rf (66) 
1 — Ao 

In other words, increasing the chain length increases the mean density uniformly, as if the floppy 
chain is simply winding around more under the same confinement. This uninteresting behavior is 
refined by two interactions that have been ignored in getting (66) . First is the stiffness of successive 
pair orientations, equivalent to the monomers having coupled orientational degrees of freedom, a 
topic that has been addressed to some extent in this format in the past and will be attended to 
more forcefully in the future. Second is the effect of non-next neighbor interactions, which we 
have here studied in a preliminary fashion. The mean field v(r) that has been enountered will also 
have the effect of correlating pair orientations, and will of course alter the nature of the long chain 
resonant state. Let us see how this works. To keep the extrapolation from (63-66) transparent, let 
us rewrite (49) (using ( = e 13 ^') as, 

Oi(C^) = Mm v ^v-r 1 v-2ln((pe-^) 1 / 2 \(p- l e^) 1 /^) 

1 



(67) 



2v -(f)- 1 -21n((pe- te ) 1 /2|^) + l n (^|/_ ( pe -bvy/2 w ( pe -/3vy/2^ 
Hence, according to (53), ip = ((~ 1 e l3v ) 1 / 2 ifj is normalized, (Y>|i/>) = 1 and from (54), 

n = Ni) 2 (68) 

Near resonance is now signaled by the approximate validity of 

(/ - (Ce-^) 1/2 u;(Ce-^) 1/2 ) ^ = A^ (69) 

with small A > 0. Eq. (50) then tells us at once that using the exact 

tjT^ = -n, (70) 
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we have 



N = l/X (71) 

There are two ways of making use of the approximation 69. Most directly, we substitute (68-71) 
into (67), obtaining 

1 



(3Q\ ~ Min n 



^n-0n-21n((Ce^ n ) 1/2 |n 1 / 2 } 



(72) 



This is obviously too sweeping an approximation: by using what is effectively a first order correction 
in the argument of a variational principle, the role of the next-neighbor w in f3Qi has vanished. 
But we can pick up the next order by working instead at the "profile equation" level. It is only 
necessary to substitute (68, 70, 71) directly into (69) to rewrite the latter as 

((e^MCe-^W = (l - 1) n 1 ^, (73) 

and hence as 

= An 1 / 2 (74) 

(fj,' is the global chemical potential). The operator in (74) (x 2 ) has a positive kernel, hence 
by Jentchke's extension of Perron-Frobenius [13], the positive eigenfunction is unique (up to a 
multiplicative constant) and A is real and maximal. This leads to what is in principle a simple 
numerical iteration: start eg. with 4> = 0, compute the eigenfunction n 1 / 2 and normalize it to 
get / n{r)dr = N . At the k th stage of the M-fold iteration, replace 4> by {k/M)<j> and repeat the 
process using the current function re. An alternative strategy is to parameterize re and determine the 
parameters by Galerkin, i.e. integrate (74) with weight function and solve the resulting algebraic 
equation. 



VI. CONCLUDING REMARKS 



In this inital study, we have made one major approximation and several simplifying restrictions. 
The approximation is of course that of mean field, or selective neglect of fluctuations. On the 
assumption that fluctuations are Gaussian to leading order (examples in which this is not the case 
are far from rare, see e.g. [14]), a correction sequence is in primciple routine: we expand InHo in 
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(37) and (42) around v of (44). Using S 2 n /6v(r)6v(r') = 6n {r)/5v(r')\ 



VI 



= exp <^ -f3 ( J) [/i - t;4>d - v] + -v ■ <f>~ 1; D + -(w - v) ■ — \v + <p~ l (v - v) + . . . 





It readily follows (see e.g. [15]) that the density profile is given to the next order by simply 
averaging the density over the Gaussian field fluctuation: 



where M(r, r') = _1 (r, r') — n' [fi — \<§>d — v — A]; here, A is the field amplitude fluctuation at r. 
The simplifying assumptions attend more directly to the physics, and these assumptions depend 
very much on the nature of the system to be studied. Taking these assumptions in order, we first 
emebedded out system in a monomer grand ensemble. Since preprocessing assures that in practice 
one does not deal with the resulting extreme polydispersity, a fixed TV ensemble is more relevant 
than fixed £. The corresponding inverse mapping has been attended to on numerous occasions (see 
e.g. [16]). The same formalism is indeed available here (hinted at in eqs.(24),(25)). 

Another restriction was to attractive long-range forces. We found however that the profile 
equation could indeed be pushed into the partially repulsive regime, although the validity of the 
minimum principle was in question; this is closely related to the functional Fourier transform for 
the repulsive component, likewise under unceirtain control. An alternative approach lies in the use 
of the mean spherical model [17, 18] and its extensions. This is the aim of ongoing research. 

Of course, there is the implicit assumption that pair forces suffice, whereas the action of pairs 
on singlets is a frequent important observation, leading e.g. to dihedral angle dependence in chains. 
Typically (see e.g. [19] for a very primitive example) one can simply create a multi-unit monomer 
to encompass only such forces, which then appear once more as pair forces. 

Most blatantly, our approach has been restricted to homopolymers. Since the set of degrees of 
freedom of a monomer can include monomer type, this is no restriction at all if one is studying 
the effect on the full population of an assumed relative frequency of occurrence of next neightbors 
hetero-pairs. However, if we attribute a sense to the chain and the AB frequency differs from the 
BA frequency, the very convenient symmetry of the operator w is lost, and with it, the possibility 
e.g. of a specific long sequence of monomer types. The case of non-symmetric w has indeed been 
studied [20] , but exercising the kind of control that we have in our current formulation remains a 




(76) 
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challenge. 
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